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Abstract. We present equation of state data of shock compressed hydrogen and deuterium. These have been calculated 
in the physical picture by using ab initio molecular dynamics simulations based on finite temperature density functional 
theory as well as in the chemical picture via the Saha-D model. The results are compared in detail with data of shock 
wave experiments obtained for condensed and gaseous precompressed hydrogen and deuterium targets in a wide range 
of shock compressions from low pressures up to megabars. 

PACS. XX.XX.XX No PACS code given 



1 Introduction 

The equation of state (EOS) of hydrogen and its isotopes has 
been in the focus of research for many years for several rea- 
sons. In models of stellar and planetary interiors |QJ|3] hydro- 
gen is the most abundant element and its EOS is the most im- 
portant component for the results. Deuterium and tritium are 
target materials in inertial confinement fusion experiments [4 |. 
Therefore, a lot of experimental and theoretical efforts were 
done to understand the behavior of hydrogen, deuterium, and 
tritium in a wide range of densities and temperatures. Recent 
developments in shock-wave experiments have enabled an ac- 
cess to a precise database in the megabar pressure range. Single 
or multiple shock-wave experiments have been performed for 
hydrogen (or deuterium) by using, e.g., high explosives []5], gas 
guns [6 1, pulsed power EH9), or high-power lasers llTOlHTI . 
The strongly correlated states of warm dense matter cover a 
wide range, from an atomic-molecular mixture at low temper- 
atures to fully ionized weakly coupled plasma at high tempera- 
tures. In particular, the characterization of the transition region 
(partially ionized plasma) is a great challenge both to theory 
and experiment since the bound states exhibit a highly tran- 
sient nature. This region of the phase diagram is, however, of 
great relevance for planetary interiors. Important experimental 
information is also gained from helioseismology data [ 12 1 that 
allows to check and correct theoretical models very accurately 
in the high temperature limit. 

Some theoretical methods yield accurate results for limit- 
ing cases which then can be used to benchmark more gen- 
eral but approximate methods. For instance, an exact asymp- 
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totic expansion of thermodynamic functions can be given in 
the limit of almost fully ionized, low density plasma 11 1 311 1 41 . 
Ab initio simulation techniques such as Path Integral Monte 
Carlo (PIMC) 115I16L quantum monte carlo (QMC) lfT7HT9l or 
finite-temperature density functional theory molecular dynam- 
ics (FT-DFT-MD) simulations B20U211 which treat quantum ef- 
fects and correlations systematically have taken a great benefit 
from the rapid progress in computing power. These methods 
provide very accurate and reliable results for a variety of prob- 
lems and systems, especially for warm dense matter. In addition 
to these approaches, advanced chemical models developed for 
partially ionized plasmas 11221 have also been applied for warm 
dense matter for a long time [23-32 1. 

In the present work we compare the results of the Saha-D 
model and of FT-DFT-MD simulations with shock-wave ex- 
periments for hydrogen and deuterium which were performed 
for different initial densities in a wide range of shock compres- 
sions. We find good agreement so that these models can also be 
used to give detailed predictions for high-pressure states that 
will be probed in future experiments by varying the initial con- 
ditions accordingly. 

Our paper is organized as follows. The FT-DFT-MD sim- 
ulations are explained in Section [2] and the Saha-D model in 
Section [3] We present our results for the Hugoniot curves in 
Section|4]and, finally, give a short summary in Section|5] 



2 FT-DFT-MD simulations 

FT-DFT-MD simulations are a powerful tool to describe warm 
dense matter [33-46 1. Correlation and quantum effects are con- 
sidered by a combination of classical molecular dynamics for 
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the ions and density functional theory for the electrons. We 
use the plane wave density functional code VASP (Vienna Ab 
Initio Simulation Package) B47H491 to perform molecular dy- 
namics simulations. VASP applies Mermin's finite tempera- 
ture density functional theory [ 50 1 which allows us to treat the 
electrons even at higher temperatures on a quantum level. Pro- 
jector augmented wave potentials [51 1 were used and we ap- 
plied a generalized gradient approximation (GGA) within the 
parameterization of PBE (52). The plane wave cutoff E cut has 
to be chosen high enough to obtain converged EOS data ll36l 
l53ll . A convergence of better than 1% is secured for E cut = 
1200 eV which was used in all calculations presented here. In 
the MD scheme of VASP the Born-Oppenheimer approxima- 
tion is used, i.e. the dynamics of the ions is treated within a 
classical MD with inter-ionic forces obtained by FT-DFT cal- 
culations via the Hellmann-Feynman theorem. The electronic 
structure calculations were performed for a static array of ions 
at each MD step. This was repeated until the EOS measures 
were converged and a thermodynamic equilibrium was reached. 

The simulations were done for 256 atoms in a supercell 
with periodic boundary conditions. A Nose thermostat [54] con- 
trolled the temperature of the ions, and the temperature of the 
electrons was fixed by Fermi weighting the occupation of the 
electronic states [48]. Sampling of the Brillouin zone using up 
to 14 k-points showed that well converged results were ob- 
tained using Baldereschi's mean value point [55] for 256 par- 
ticles. The same convergence behavior has previously been re- 
ported for water [56]. The size of the simulated supercell fixed 
the density of the system. The internal energy was corrected 
due to zero point vibrations of the molecules at low tempera- 
tures, taking into account quantum contributions of a harmonic 
oscillator for each molecule [57]. For this procedure the num- 
ber of molecules has to be known, which was obtained by eval- 
uating the pair correlation function, see Ref. Il53l for details. 
The system was simulated 1000-1500 steps further after reach- 
ing the thermodynamic equilibrium to ensure a small statistical 
error. The EOS data were then obtained by averaging over all 
particles and simulation steps in equilibrium. 



3 Saha-D model 

The Saha-D model EOS is based on the chemical picture El 
30 58 1 which represents the plasma as a mixture of interacting 
electrons, ions, atoms, and molecules. We consider the follow- 
ing components for hydrogen and deuterium: e~, A, A + , A2, A£ , 
(A : H, D). For this case the Helmholtz free energy reads: 



Table 1. Parameters of A2 - A2, At - A£, A - A repulsion in 
Eq. 0; flo is the Bohr radius. 



F({Nj}, V, T) = J] Ff> + Ff> + AFf' } + AFf f) . 



(1) 



We shortly outline the approximations in which these three 
contributions were treated; for details, see [30.59]. The first two 
terms of Eq. (0 are the ideal gas contributions of heavy parti- 
cles (atoms, ions, and molecules) and of electrons. The latter 
corresponds to the partially degenerate ideal Fermi gas [60 1. 
The last two terms describe corrections due to Coulomb in- 
teractions and short-range interactions between heavy parti- 
cles. The Coulomb interaction effects of charged particles are 
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considered here within a modified pseudopotential approach 
EBEDI62. 

The electron-electron and ion-ion interaction is each treated 
by using the Coulomb potential. For the effective electron-ion 
interaction we apply a pseudopotential using the Glauberman 
form [63 1 . The parameters of the pseudopotentials and of the 
corresponding pair correlation functions were determined from 
the general conditions of local electroneutrality and dipole screen- 
ing, from the non-negativity constraint for the pair correlation 
functions, and from a relation between the screening cloud am- 
plitude and the depth of the electron-ion pseudopotential. In the 
weak coupling limit, this approximation coincides with the De- 
bye model but in the strong coupling limit it is much softer and 
demonstrates a quasi-crystalline behavior. 

At high densities as typical for shock-compressed hydro- 
gen the short-range repulsion between composite heavy parti- 
cles (A, A2, At) becomes very important. This effect is taken 
into account in the Saha-D model within a simple soft-sphere 
approximation [64| which is modified for a mixture of soft 
spheres with different radii. In this case the effective packing 
fraction Y is calculated via the individual diameters cr ; of each 
particle species in correspondence with the one-fluid approxi- 
mation: 



nn<j r 



p3\ 



1/3 



(2) 



(jj is the diameter of the soft spheres in the respective poten- 
tial Vss(r) = e(r/crj)~ s . The contribution of the intermolecu- 
lar repulsion dominates the EOS of dense hydrogen and deu- 
terium in a wide range of pressures at low temperatures. The 
contribution of atom-atom and atom-molecule repulsion be- 
comes important at elevated temperatures. The parameters for 
the soft-sphere repulsion for A2 - A2, A2 - A, and A - A are 
chosen according to the spherically symmetric parts of the ef- 
fective interaction potentials of the non-empirical atom- atom 
approximation ll65l . The key parameter of this approximation is 
the ratio of corresponding soft-sphere diameters for atoms and 
molecules, o~a/(Ta 2 - This ratio determines the change of intrin- 
sic volumes of two atoms in comparison with that of a molecule 
(2Va/Va 2 )', it determines the effective shift of the dissociation 
and ionization equilibrium in warm dense hydrogen and deu- 
terium. All parameters of the soft-sphere repulsion are given in 
Table The parameter e is chosen such that our soft-sphere 
potential for molecule-molecule repulsion will be close to the 
potential 1631 at a distance r = 2a (in this case e = 0.138eV). 

In order to take into account the existence of condensed 
states (liquid and solid) for hydrogen and deuterium, the at- 
traction term in the free energy has to be considered together 
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with the soft-sphere repulsion. In our case AF^ n ' ) reads: 
AF ( '" r) = AF ss +AF Mr , 



(3) 



AF„ 



The attractive corrections [66] are independent of temperature. 
The choice of 6 = 1 as in our case corresponds to a van der 
Waals-like approximation. The parameter B supplies the cor- 
rect sublimation energy of a molecular system in the condensed 
state. 



4 Results 

We calculated Hugoniot curves based on the FT-DFT-MD and 
Saha-D EOS data sets for hydrogen and deuterium. In Fig. [T]we 
compare our calculations with gas-gun experiments [67. 68] on 
liquid hydrogen and deuterium. The experiments for deuterium 
and hydrogen are reproduced by the Saha-D model within the 
uncertainties of the experiments. The FT-DFT-MD Hugoniot 
curves reproduce the experiments with less precision. In the 
case of hydrogen the compression rate is slightly overestimated 
above 5 GPa. For deuterium this occurs for pressures above 
15 GPa. Furthermore, the FT-DFT-MD curve lies above the ex- 
periments for pressures below 10 GPa. Between 10 GPa and 
20 GPa the experimental data are reproduced within error bars. 




Fig. 1. Shock compression of liquid hydrogen (blue) and deu- 
terium (red). The Hugoniot curves obtained with the Saha-D 
model (solid line) and the FT-DFT-MD (dashed line) are com- 
pared with shock wave experiments of Nellis et al. [67| (cir- 
cles) and Holmes et al. [68| (squares) 



The less precision of the FT-DFT-MD data in the region 
of the gas gun experiments is connected with the relatively 
abrupt onset of dissociation processes which lead to an increase 
of the compression ratio at about 20 GPa and 4200 K as it 
has already been reported ll69l . The reason for this behavior 
can be related to the underestimation of the fundamental band 
gap in DFT-based electronic structure calculations. More ac- 
curate exchange-correlation functionals than PBE, specifically 
derived for finite temperatures, are urgently needed for warm 
dense matter. QMC calculations treat the exchange-correlation 
directly and are not affected by this approximation. Recent QMC 
calculations find in fact a shift of the dissociation region com- 
pared to DFT [19 1 but show no different results for the EOS 



for conditions relevant for planetary interiors [18. 70 . For in- 
stance, the Saha-D model yields 10% dissociation at 60 GPa 
and a temperature of 13000 K along the deuterium Hugoniot 
curve. Along the hydrogen Hugoniot curve dissociation occurs 
at about 15 GPa and 3000 K within the FT-DFT-MD model. 
The Saha-D model predicts 10% dissociation at 50 GPa and a 
temperature of 14000 K along the hydrogen Hugoniot curve. 

The presented theoretical predictions both agree with the 
experiments within 10% accuracy in compression ratio and can 
therefore describe the principal behavior of the obtained results 
well. The compression reached in shock compressed hydrogen 
is higher than in deuterium at the same pressure. This is repro- 
duced by the calculated Hugoniot curves, see also [45|. 

The FT-DFT-MD hydrogen EOS data were also used to cal- 
culate the deuterium Hugoniot curve. To adjust for the deu- 
terium initial conditions we considered the initial conditions of 
the hydrogen EOS at half the deuterium density given in the 
experiment (0.171 g/cm 3 ), i.e. 0.0855 g/cm 3 . Plotting the re- 
sulting pressure versus the compression ratio, both Hugoniot 
curves are almost the same. On the other hand, calculating the 
pressure with the deuterium EOS and adjusting for the hydro- 
gen initial conditions in the same way, the resulting Hugoniot 
curves are identical within a smaller error than the statistical er- 
ror of the FT-DFT-MD simulations. The different compression 
ratios of hydrogen and deuterium as seen in the experiment are, 
therefore, only slightly caused by differences in the EOS data 
of both isotopes at warm dense matter conditions. The differ- 
ence in the compression ratios is mainly due to the fact that the 
densities of the liquid targets at 20 K do not scale exactly by a 
factor of two. Scaling the deuterium density to that of hydrogen 
the initial density of liquid deuterium would be 0.0855 g/cm 3 
which differs from the value relevant for liquid hydrogen which 
is 0.071 g/cm 3 . This deviation of about 20% entails the differ- 
ent Hugoniot curves. 
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Fig. 2. Shock compression of liquid hydrogen (blue) and deu- 
terium (red). Temperatures along the Hugoniot curves as ob- 
tained with the Saha-D model (solid line) and FT-DTF-MD 
(dashed line) are displayed as function of pressure and com- 
pared with shock-wave experiments (squares) 



The temperature along the Hugoniot curves as predicted by 
the theoretical models is shown in Fig. [2] and compared with 
gas-gun experimental data [68 1. The temperatures measured 
in the experiments are in general higher than predicted by the 
Saha-D model and the FT-DFT-MD simulations, except for two 
deuterium points above 22 GPa which are below the Saha-D 
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curve. Again, onset of dissociation causes the slight kink in the 
FT-DFT-MD curves at about 18 GPa (D 2 ) and 15 GPa (H 2 ). 
The general behavior indicated by the experiments can be re- 
produced: the temperature along the hydrogen Hugoniot curve 
is higher than that for deuterium at the same pressure. The max- 
imum deviation of both theoretical models from the experimen- 
tal data is about 400 K. 

We have also applied both theoretical EOS data sets to cal- 
culate the Hugoniot curves for different initial conditions in 
order to study the compression behavior of the hydrogen iso- 
topes for a wide range of densities off the principal Hugoniot. 
Fig. [3] shows Hugoniot curves with respect to initial condi- 
tions as chosen in recent shock-wave experiments with precom- 
pressed targets 07111721 . Two experiments were performed with 
gaseous targets at 1.5 kbar (po = 0.1335 g/cm 3 ) and 2.0 kbar 
(po = 0.153 g/cm 3 ). Two other data sets were obtained with liq- 
uid (po = 0.171 g/cm 3 ) and solid (po - 0.199 g/cm 3 ) deuterium 
targets. 
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Fig. 3. Deuterium single shock principal Hugoniot curves 
starting from different initial densities (gaseous, liquid, and 
solid) as derived from FT-DFT-MD (black), Saha-D (35) (red), 
RPIMC G2 (blue), DPIMC d (dark red). Dotted line: p = 
0.1335 g/cm 3 , dashed line: po = 0.153 g/cm 3 , solid line: po = 
0.171 g/cm 3 , dot-dashed line: po = 0.199 g/cm 3 . Shock wave 
experiments: Nellis et al. |67| (blue circles), Holmes et al. |68| 
(blue squares), Knudson et al. [7, 8| (open green circles), Gr- 
ishechkin et al. iPTD (red triangles), Boriskov et al. [72 1 (red 
squares and diamonds). The arrows at the top show the limiting 
compression for ultra-high pressures (4 x po) for each principal 
Hugoniot. 

The two theoretical results and the experiments show the 
same general behavior: the attained absolute density is higher 
the more precompressed the target is. It has to be pointed out 
that the maximum compression ratio shows the inverse behav- 
ior, it decreases with higher precompression. Even so, the max- 
imum density that can be probed in single shock experiments 
increases with precompression. The theoretical predictions of 
the two methods for the maximum density agree at the con- 
ditions of the experiments with gaseous and liquid targets and 
range from 0.65 g/cm 3 to 0.775 g/cm 3 . For the initial condition 



in the solid, there is a slight difference: FT-DFT-MD predicts 
0.83 g/cm 3 and Saha-D 0.87 g/cm 3 . The pressure at this maxi- 
mum compression density is also different for the two models; 
it ranges from 30 to 50 GPa within FT-DFT-MD and from 80 
to 150 GPa according to Saha-D. These values cannot be dis- 
criminated via the few experimental points. 

Fig.@]shows the temperatures along the Hugoniot curves of 
deuterium using the initial conditions of the experiments with 
liquid and precompressed gaseous targets [68. 7111751. 
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Fig. 4. Temperature of shock compressed deuterium for differ- 
ent initial conditions predicted by FT-DFT-MD (black), Saha- 
D 132) (red), RPIMC (73) (blue), DPIMC (71 (dark red), 
and the asymptotically strict Saha-S [59 1 (orange) model in 
comparison with experiments using liquid [Holmes et al. [68 1 
(blue squares) and Bailey et al. 11751 (green hexagons)] and 
gaseous targets [Grishechkin et al. fT\\ (red stars)]. Dotted, 
dashed, and solid lines correspond to initial deuterium densi- 
ties of po = 0.1335,0.153,0.171 g/cm 3 , respectively. Arrows 
indicate the limiting compression for ultra-high pressures for 
each Hugoniot curve. 

We have to note again that with a higher precompression 
also a higher density can be reached. The measurements in- 
dicate a temperature of about 2 eV for all three initial condi- 
tions. These results can be reproduced by the theoretical mod- 
els within the error bars. The densest states are reached for the 
liquid deuterium target. The experiments of Bailey et al. f75l 
show a maximum density between 0.7 and 0.8 g/cm 3 which is 
reproduced by both theories. The temperatures measured in the 
experiments are underestimated by FT-DFT-MD and overesti- 
mated by Saha-D, with increasing deviation for higher temper- 
atures. The limit of the Saha-D model at high temperatures can 
be checked by comparing with results obtained by the Saha-S 
model [59] which is asymptotically exact for To <k l,«/t 3 «: 1. 
Such a comparison shows that Saha-D (together with PIMC [ 73 
1741 ) yields the correct high-temperature limit. The deviation of 
the Saha-S from the Saha-D curves at lower temperatures is 
due to the fact that the Saha-S model is no longer applicable 
for these parameters. In particular, the Saha-S model does not 
take into account the short-range repulsion effects of compos- 
ite particles (A,A2,AJ). Interestingly, the temperature at max- 
imum compression is almost independent of the initial condi- 
tions. The predictions of the theoretical models show that the 
curves are shifted only to higher densities while the tempera- 
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ture remains almost constant. This finding is supported by the 
experimental results. 

5 Conclusions 

We have calculated the EOS of deuterium and hydrogen with 
FT-DFT-MD simulations in the physical picture and the Saha- 
D model in the chemical picture over a wide range of density 
and temperature which enabled us to compare those results 
with recent shock-wave experiments. The theories predict, in 
agreement with the experiments, that higher compressed den- 
sities can be reached using precompressed targets, while the 
maximum compression ratio decreases. We compare also the 
temperature along the Hugoniot curve of deuterium with exper- 
imental data and find that only the density is affected by pre- 
compression, while the temperature remains almost the same 
along the different pathways. This leads to an increased pres- 
sure with higher precompression along the Hugoniot. 

A check of EOS models against experiments within the 
WDM regime is available mostly for the relatively limited den- 
sity and temperature range along the principal Hugoniot curve. 
Experiments producing shock waves within precompressed tar- 
gets enable to check the quality of EOS models in a wider range 
in the phase diagram. Both EOS (based on FT-DFT-MD simu- 
lations and the Saha-D model) could reproduce the experimen- 
tal data. On the other hand, neither the Saha-D model, which 
uses effective two-particle potentials with parameters that have 
been chosen to match physical constraints, nor the FT-DFT- 
MD method, which has no adjustable parameters, can repro- 
duce all experimental features precisely. Nevertheless, exper- 
imental data is still not available in the needed quantity and 
precision to allow for a definite decision of which model has 
to be used to describe all quantities in all ranges of the phase 
diagram. We therefore look forward to future high-pressure ex- 
periments, especially off the principal Hugoniot curve. 

The combination of an advanced chemical model with an 
ab initio approach yields a reasonable description of warm dense 
hydrogen because the low-density molecular liquid, the strongly 
correlated warm dense fluid, as well as the hot plasma can be 
described adequately, see also Il76l . Simultaneously this com- 
bination saves computational power as the treatment of a low- 
density molecular liquid is increasingly demanding when using 
FT-DFT-MD simulations. The treatment of a free energy min- 
imiztion model like Saha-D is much less expensive regarding 
computational time. Accordingly, this combination provides an 
opportunity to construct a wide-range EOS for planetary inte- 
rior modelling for which a database from ambient conditions 
up to pressures of several tens of megabar and temperatures 
up to about 100.000 K is needed. This project is going to be 
compiled for and will be be applied to model the interior of 
Jupiter 11771 . 
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